Ejercicio metodológico de interpretación de puntos en el marco de SIMOCUTE

Resumen

FIESTA es una biblioteca del lenguaje de programación R para análisis de datos de inventarios forestales basados en muestras. En este documento se explica el uso de FIESTA mediante un caso de ejemplo basado en un conjunto de datos de parcelas de monitoreo, ubicadas en la zona norte de Costa Rica.

Introducción

FIESTA (Forest Inventory Estimation and Analysis) es una biblioteca de software, desarrollada en el lenguaje de programación R, para el análisis de datos de inventarios forestales basados en muestras. Es desarrollada por el Programa de Inventario y Análisis Forestal (FIA, por sus siglas en inglés) del Servicio Forestal del Departamento de Agricultura de los Estados Unidos.

Este documento describe el procesamiento, mediante FIESTA, de un conjunto de datos de 26125 puntos de muestreo de 1045 parcelas de monitoreo, ubicadas en la zona norte de Costa Rica.

Recursos

• Forest Inventory Estimation and Analysis (FIESTA) (sitio oficial de FIESTA)
• FIESTA - Photo-Based Module (módulo de inventarios basados en fotografías de FIESTA)

Trabajo previo

Instalación de software

Sistema base de R y herramientas de desarrollo

Para trabajar con FIESTA, se debe instalar:

• El sistema base de R. R es un lenguaje de programación enfocado en análisis estadístico y visualización de datos.
• La interfaz de desarrollo integrada RStudio Desktop, la cual proporciona un editor de texto y otras herramientas para escribir programas en R y visualizar sus resultados, entre otras facilidades.

Versión mínima de R

De acuerdo con su documentación, la versión mínima de R que requiere FIESTA, a la fecha de escritura de este documento (2023-12-02), es la 4.2.0. Puede consultar la versión que tiene instalada al ejecutar el siguiente comando en la consola de R:

# Versión de R
R.version.string

[1] "R version 4.3.2 (2023-10-31 ucrt)"

La salida del comando anterior debe indicar que su versión de R es mayor o igual a 4.2.0.

Paquetes de R

FIESTA

El paquete FIESTA está disponible en CRAN (Comprehensive R Archive Network), un repositorio en línea que alberga una amplia colección de paquetes y extensiones para el lenguaje de programación R, lo que facilita su instalación y actualización a nuevas versiones.

Para usar FIESTA, debe instalarse primero en una interfaz de R. Puede utilizar la función install.packages().

# Instalación del paquete FIESTA
install.packages("FIESTA")

El resultado de la instalación puede verificarse al cargar el paquete con la función library().

# Carga del paquete FIESTA
library(FIESTA)

Si el comando anterior no genera ningún mensaje de error, FIESTA debe haberse instalado adecuadamente.

Otros

Además de FIESTA, se recomienda instalar los siguientes paquetes para procesamiento de datos y visualización de resultados.

# Paquete para el desarrollo de documentos computacionales
install.packages("rmarkdown")

# Colección de paquetes para análisis de datos
install.packages("tidyverse")

# Estilos para gráficos de tidyverse
install.packages("ggthemes")

# Paquete para limpieza de datos
install.packages("janitor")

# Paquete para tablas interactivas
install.packages("DT")

# Paquete para graficación interactiva
install.packages("plotly")

# Paquete para mapas interactivos
install.packages("leaflet")

# Funciones adicionales para leaflet
install.packages("leaflet.extras")

# Funciones adicionales para leaflet
install.packages("leafem")

Luego de instalarlos, debe cargar los paquetes con la función library().

# Carga de paquetes adicionales
library(rmarkdown)
library(tidyverse)
library(ggthemes)
library(janitor)
library(DT)
library(plotly)
library(leaflet)
library(leaflet.extras)
library(leafem)
library(sf) # se instala con FIESTA

Obtención de este repositorio

Este documento forma parte de un repositorio en GitHub, una plataforma en línea para compartir código fuente de aplicaciones, basada en el sistema de control de versiones Git. El repositorio contiene el código fuente del documento y los datos que se utilizan en los ejemplos. Su dirección es https://github.com/mesa-monitoreo-puntos/fiesta.

Puede descargar el repositorio, como un archivo ZIP, de https://github.com/mesa-monitoreo-puntos/fiesta/archive/refs/heads/main.zip

También puede “clonarlo” mediante el comando clone de Git:

# Clonación de este repositorio
git clone https://github.com/mesa-monitoreo-puntos/fiesta.git

Una vez que el repositorio haya sido descargado o clonado, puede abrirse con RStudio o con otra herramienta de desarrollo.

Variables generales

En esta sección se definen algunas variables generales del proceso.

Código para la definición de variables generales

# Ruta a los datos de puntos de muestreo
ARCHIVO_PUNTOS <- "datos/Resul_Fin_2Grupos_SinReplicas_csv_Fix.csv"

# Colores
COLOR_VEGETACION         = "#00a600" # Corine EU - Coniferous forest
COLOR_SIN_VEGETACION     = "#dcdcc8" # Corine CR - Tierras desnudas o degradadas
COLOR_AGUA               = "#00ccf2" # Corine EU - Water courses
COLOR_NUBES_Y_SOMBRAS    = "#800080" # Púrpura
COLOR_SIN_INFORMACION    = "#000000" # Negro

COLOR_MC_BOSQUES         = "#266900" # Corine CR - Bosque denso
COLOR_AGRICULTURA        = "#Becd05" # Corine CR - Mosaico de cultivos
COLOR_GANADERIA_Y_PASTOS = "#ffffa6" # Corine CR - Pastos limpios
COLOR_ZONAS_HUMEDALES    = "#00ccf2" # Corine EU - Water courses
COLOR_NO_CLASIFICACION   = "#000000" # Negro
COLOR_INFRAESTRUCTURA    = "#f6d9df" # Corine CR - Zona urbana continua
COLOR_OTRAS_TIERRAS      = "#005acf" # Corine CR - Canales

COLOR_NULO               = "#808080" # Gris

# Paleta de colores de tipos de cobertura de la tierra
COLORES_COBERTURA <- 
    c(
        "1000-Vegetacion"         = COLOR_VEGETACION,
        "T1-1000-Vegetacion"      = COLOR_VEGETACION,
        "T2-1000-Vegetacion"      = COLOR_VEGETACION,
        "2000-Sin vegetacion"     = COLOR_SIN_VEGETACION, 
        "T1-2000-Sin vegetacion"  = COLOR_SIN_VEGETACION,
        "T2-2000-Sin vegetacion"  = COLOR_SIN_VEGETACION,
        "3000-Agua"               = COLOR_AGUA, 
        "T1-3000-Agua"            = COLOR_AGUA,
        "T2-3000-Agua"            = COLOR_AGUA,
        "4000-Nubes y sombras"    = COLOR_NUBES_Y_SOMBRAS,
        "T1-4000-Nubes y sombras" = COLOR_NUBES_Y_SOMBRAS,
        "T2-4000-Nubes y sombras" = COLOR_NUBES_Y_SOMBRAS,
        "5000-Sin informacion"    = COLOR_SIN_INFORMACION,
        "T1-5000-Sin informacion" = COLOR_SIN_INFORMACION,
        "T2-5000-Sin informacion" = COLOR_SIN_INFORMACION,
        "Nulo"                    = COLOR_NULO,
        "T1-Nulo"                 = COLOR_NULO,
        "T2-Nulo"                 = COLOR_NULO
    )

# Paleta de colores de tipos de uso de la tierra
COLORES_USO <- 
    c(
        "1000-MC bosques"       = COLOR_MC_BOSQUES,
        "T1-1000-MC bosques"    = COLOR_MC_BOSQUES,
        "T2-1000-MC bosques"    = COLOR_MC_BOSQUES,
        "2000-Agricultura"      = COLOR_AGRICULTURA,
        "T1-2000-Agricultura"   = COLOR_AGRICULTURA,
        "T2-2000-Agricultura"   = COLOR_AGRICULTURA,
        "3000-Ganad y past"     = COLOR_GANADERIA_Y_PASTOS,
        "T1-3000-Ganad y past"  = COLOR_GANADERIA_Y_PASTOS,
        "T2-3000-Ganad y past"  = COLOR_GANADERIA_Y_PASTOS,
        "4000-Zonas humed"      = COLOR_ZONAS_HUMEDALES,
        "T1-4000-Zonas humed"   = COLOR_ZONAS_HUMEDALES,
        "T2-4000-Zonas humed"   = COLOR_ZONAS_HUMEDALES,
        "5000-Infraest"         = COLOR_INFRAESTRUCTURA,
        "T1-5000-Infraest"      = COLOR_INFRAESTRUCTURA,
        "T2-5000-Infraest"      = COLOR_INFRAESTRUCTURA,
        "6000-Otras tierras"    = COLOR_OTRAS_TIERRAS,
        "T1-6000-Otras tierras" = COLOR_OTRAS_TIERRAS,
        "T2-6000-Otras tierras" = COLOR_OTRAS_TIERRAS,
        "7000 No clasif"        = COLOR_NO_CLASIFICACION,
        "T1-7000 No clasif"     = COLOR_NO_CLASIFICACION,
        "T2-7000 No clasif"     = COLOR_NO_CLASIFICACION,
        "Nulo"                  = COLOR_NULO,
        "T1-Nulo"               = COLOR_NULO,
        "T2-Nulo"               = COLOR_NULO
    )

# Área de estudio en hectáreas
AREA_ESTUDIO = 318167.67

Carga y limpieza de datos

El conjunto de datos que se utiliza en este documento para ilustrar el uso de FIESTA es el producto de la interpretación de 1045 parcelas en la zona norte de Costa Rica. Cada parcela contiene 25 puntos, por lo que el total de puntos es 26125. Cada punto se interpretó, mediante imágenes satelitales, en dos tiempos:

• Tiempo 1 (t1): entre 2005 y 2007.
• Tiempo 2 (t2): en 2019.

La interpretación fue realizada por un equipo de 26 personas de diferentes instituciones y organizaciones.

Los datos de puntos de muestreo se cargan de un archivo CSV y los nombres de las columnas se “limpian” para evitar números al inicio y otros problemas que dificultan su manejo.

Código para la carga y limpieza de datos

# Carga de datos de puntos de muestreo
puntos <- read_delim(ARCHIVO_PUNTOS)

# Limpieza de los nombres de columnas
puntos <- clean_names(puntos)

# Reemplazo de valores NA por la hilera "Nulo"
puntos <- 
    puntos |> 
    mutate(
        t1_cobertura = replace_na(t1_cobertura, "Nulo"),
        t2_cobertura = replace_na(t2_cobertura, "Nulo"),
        t1_uso = replace_na(t1_uso, "Nulo"),
        t1_uso = replace_na(t1_uso, "Nulo")
    )

# Coversión a de t1_cobertura, t2_cobertura, t1_uso y t2_uso a factores
puntos <-
    puntos |>
    mutate(
        t1_cobertura = factor(t1_cobertura, levels = unique(t1_cobertura)),
        t2_cobertura = factor(t2_cobertura, levels = unique(t2_cobertura)),
        t1_uso = factor(t1_uso, levels = unique(t1_uso)),
        t2_uso = factor(t2_uso, levels = unique(t2_uso))        
    )

Código para el despliegue de la tabla

# Despliegue de los datos de puntos de muestreo en una tabla
puntos |>
    select(plot_id, sample_id, t1_cobertura, t2_cobertura, t1_uso, t2_uso) |>
    datatable(
        caption = "Puntos de muestreo de cobertura y uso de la tierra",
        rownames = FALSE,
        colnames = c(
            "plot_id", "sample_id", 
            "t1_cobertura", "t2_cobertura",
            "t1_uso", "t2_uso"
        ),
        options = list(
            pageLength = 5,
            language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
        )
    )

Análisis con paquetes clásicos de R

Cobertura de la tierra

Mapa

Código para el despliegue del mapa

# Conversión de datos de puntos de muestreo a objeto sf (vectorial de puntos)
geo_puntos <-
    puntos |>
    select(plot_id, sample_id, lon, lat, t1_cobertura, t2_cobertura) |>
    st_as_sf(
        coords = c("lon", "lat"),
        crs = 4326
  )

colores <- 
    colorFactor(
        palette = c(
            COLOR_VEGETACION, COLOR_SIN_VEGETACION, COLOR_AGUA, 
            COLOR_NUBES_Y_SOMBRAS, COLOR_SIN_INFORMACION, COLOR_NULO
        ), 
    levels = c(
        "1000-Vegetacion", "2000-Sin vegetacion", "3000-Agua",
        "4000-Nubes y sombras", "5000-Sin informacion", "Nulo"
    )
)

leaflet() |>
    addTiles(group = "OSM") |>
    addProviderTiles(
        provider = providers$Esri.WorldImagery, 
        group = "ESRI World Imagery"
    ) |>
    addProviderTiles(
        provider = providers$CartoDB.DarkMatter,
        group = "Dark Matter"
    ) |>        
    addCircleMarkers(
        data = geo_puntos,
        radius = 2,
        fillColor = ~colores(geo_puntos$t1_cobertura),
        color = ~colores(geo_puntos$t1_cobertura),      
        # clusterOptions = markerClusterOptions(),
    popup = paste(
      paste0("<strong>Parcela: </strong>", geo_puntos$plot_id),
      paste0("<strong>Muestra: </strong>", geo_puntos$sample_id),
      paste0("<strong>Cobertura en t1: </strong>", geo_puntos$t1_cobertura),
      paste0("<strong>Cobertura en t2: </strong>", geo_puntos$t2_cobertura),
      sep = '<br/>'
    ),      
        group = "Cobertura en t1"
    ) |>
  addLegend(
    position = "bottomleft",    
    pal = colores,
    values = geo_puntos$t1_cobertura,
    title = "Cobertura en t1",
    group = "Cobertura en t1"    
  ) |>      
    addCircleMarkers(
        data = geo_puntos,
        radius = 2,
        fillColor = ~colores(geo_puntos$t2_cobertura),
        color = ~colores(geo_puntos$t2_cobertura),      
        # clusterOptions = markerClusterOptions(),
    popup = paste(
      paste0("<strong>Parcela: </strong>", geo_puntos$plot_id),
      paste0("<strong>Muestra: </strong>", geo_puntos$sample_id),
      paste0("<strong>Cobertura en t1: </strong>", geo_puntos$t1_cobertura),
      paste0("<strong>Cobertura en t2: </strong>", geo_puntos$t2_cobertura),
      sep = '<br/>'
    ),      
        group = "Cobertura en t2"
    ) |>    
  addLegend(
    position = "bottomleft",    
    pal = colores,
    values = geo_puntos$t2_cobertura,
    title = "Cobertura en t2",
    group = "Cobertura en t2"    
  ) |>  
    addLayersControl(
        baseGroups = c("OSM", "ESRI World Imagery", "Dark Matter"),
        overlayGroups = c("Cobertura en t1", "Cobertura en t2"),
    )

Distribución de puntos de muestreo en t1

Código para el despliegue del gráfico de barras

# Gráfico de barras ggplot2
grafico_barras <-
    puntos |>
    ggplot(aes(x = fct_infreq(t1_cobertura), fill = fct_infreq(t1_cobertura))) +
    geom_bar(
    aes(
      text = paste0(
        "Cantidad de puntos: ", after_stat(count)
      )     
    )
    ) +
    scale_fill_manual(values = COLORES_COBERTURA) +
    xlab("Cobertura") +
    ylab("Cantidad de puntos") +
    labs(fill = "Cobertura") +
    theme_clean() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.title = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank()
  )

# Gráfico de barras plotly
ggplotly(grafico_barras, tooltip = "text") |> 
  config(locale = 'es') |>
  layout(
    showlegend = TRUE,
    legend = list(
      title = list(text = ''),
      bordercolor = 'rgba(0,0,0,0)' # Establece el color del borde a transparente
    )
  ) 

Código para el despliegue del gráfico de pastel

# Total de puntos
conteo <- puntos |>
    count(t1_cobertura)

# Asegurarse de que los colores se asignen en el orden correcto
colores_para_plotly <- COLORES_COBERTURA[as.character(conteo$t1_cobertura)]

# Gráfico de pastel plotly
plot_ly(
    conteo, 
    labels = ~ t1_cobertura, 
    values = ~ n,
    marker = list(colors = colores_para_plotly)
) |>
add_pie()

Distribución de puntos de muestreo en t2

Código para el despliegue del gráfico de barras

# Gráfico de barras ggplot2
grafico_barras <-
    puntos |>
    ggplot(aes(x = fct_infreq(t2_cobertura), fill = fct_infreq(t2_cobertura))) +
    geom_bar(
    aes(
      text = paste0(
        "Cantidad de puntos: ", after_stat(count)
      )     
    )
    ) +
    scale_fill_manual(values = COLORES_COBERTURA) +
    xlab("Cobertura") +
    ylab("Cantidad de puntos") +
    labs(fill = "Cobertura") +
    theme_clean() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.title = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank()
  )

# Gráfico de barras plotly
ggplotly(grafico_barras, tooltip = "text") |> 
  config(locale = 'es') |>
  layout(
    showlegend = TRUE,
    legend = list(
      title = list(text = ''),
      bordercolor = 'rgba(0,0,0,0)' # Establece el color del borde a transparente
    )
  ) 

Código para el despliegue del gráfico de pastel

# Total de puntos
conteo <- puntos |>
    count(t2_cobertura)

# Asegurarse de que los colores se asignen en el orden correcto
colores_para_plotly <- COLORES_COBERTURA[as.character(conteo$t2_cobertura)]

# Gráfico de pastel plotly
plot_ly(
    conteo, 
    labels = ~ t2_cobertura, 
    values = ~ n,
    marker = list(colors = colores_para_plotly)
) |>
add_pie()

Uso del suelo

Mapa

Código para el despliegue del mapa

# Conversión de datos de puntos de muestreo a objeto sf (vectorial de puntos)
geo_puntos <-
    puntos |>
    select(plot_id, sample_id, lon, lat, t1_uso, t2_uso) |>
    st_as_sf(
        coords = c("lon", "lat"),
        crs = 4326
  )

colores <- 
    colorFactor(
        palette = c(
            COLOR_MC_BOSQUES, COLOR_AGRICULTURA, COLOR_GANADERIA_Y_PASTOS, 
            COLOR_ZONAS_HUMEDALES, COLOR_INFRAESTRUCTURA, COLOR_OTRAS_TIERRAS, 
            COLOR_NO_CLASIFICACION, COLOR_NULO
        ), 
    levels = c(
        "1000-MC bosques", "2000-Agricultura", "3000-Ganad y past",
        "4000-Zonas humed", "5000-Infraest", "6000-Otras tierras",
        "7000 No clasif", "Nulo"
    )
)

leaflet() |>
    addTiles(group = "OSM") |>
    addProviderTiles(
        provider = providers$Esri.WorldImagery, 
        group = "ESRI World Imagery"
    ) |>
    addProviderTiles(
        provider = providers$CartoDB.DarkMatter,
        group = "Dark Matter"
    ) |>        
    addCircleMarkers(
        data = geo_puntos,
        radius = 2,
        fillColor = ~colores(geo_puntos$t1_uso),
        color = ~colores(geo_puntos$t1_uso),        
        # clusterOptions = markerClusterOptions(),
    popup = paste(
      paste0("<strong>Parcela: </strong>", geo_puntos$plot_id),
      paste0("<strong>Muestra: </strong>", geo_puntos$sample_id),
      paste0("<strong>Uso en t1: </strong>", geo_puntos$t1_uso),
      paste0("<strong>Uso en t2: </strong>", geo_puntos$t2_uso),
      sep = '<br/>'
    ),      
        group = "Uso en t1"
    ) |>
  addLegend(
    position = "bottomleft",    
    pal = colores,
    values = geo_puntos$t1_uso,
    title = "Uso en t1",
    group = "Uso en t1"    
  ) |>      
    addCircleMarkers(
        data = geo_puntos,
        radius = 2,
        fillColor = ~colores(geo_puntos$t2_uso),
        color = ~colores(geo_puntos$t2_uso),        
        # clusterOptions = markerClusterOptions(),
    popup = paste(
      paste0("<strong>Parcela: </strong>", geo_puntos$plot_id),
      paste0("<strong>Muestra: </strong>", geo_puntos$sample_id),
      paste0("<strong>Uso en t1: </strong>", geo_puntos$t1_uso),
      paste0("<strong>Uso en t2: </strong>", geo_puntos$t2_uso),
      sep = '<br/>'
    ),      
        group = "Uso en t2"
    ) |>    
  addLegend(
    position = "bottomleft",    
    pal = colores,
    values = geo_puntos$t2_uso,
    title = "Uso en t2",
    group = "Uso en t2"    
  ) |>  
    addLayersControl(
        baseGroups = c("OSM", "ESRI World Imagery", "Dark Matter"),
        overlayGroups = c("Uso en t1", "Uso en t2"),
    )

Distribución de puntos de muestreo en t1

Código para el despliegue del gráfico de barras

# Gráfico de barras ggplot2
grafico_barras <-
    puntos |>
    ggplot(aes(x = fct_infreq(t1_uso), fill = fct_infreq(t1_uso))) +
    geom_bar(
    aes(
      text = paste0(
        "Cantidad de puntos: ", after_stat(count)
      )     
    )
    ) +
    scale_fill_manual(values = COLORES_USO) +
    xlab("Uso") +
    ylab("Cantidad de puntos") +
    labs(fill = "Uso") +
    theme_clean() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.title = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank()
  )

# Gráfico de barras plotly
ggplotly(grafico_barras, tooltip = "text") |> 
  config(locale = 'es') |>
  layout(
    showlegend = TRUE,
    legend = list(
      title = list(text = ''),
      bordercolor = 'rgba(0,0,0,0)' # Establece el color del borde a transparente
    )
  ) 

Código para el despliegue del gráfico de pastel

# Total de puntos
conteo <- puntos |>
    count(t1_uso)

# Asegurarse de que los colores se asignen en el orden correcto
colores_para_plotly <- COLORES_USO[as.character(conteo$t1_uso)]

# Gráfico de pastel plotly
plot_ly(
    conteo, 
    labels = ~ t1_uso, 
    values = ~ n,
    marker = list(colors = colores_para_plotly)
) |>
add_pie()

Distribución de puntos de muestreo en t2

Código para el despliegue del gráfico de barras

# Gráfico de barras ggplot2
grafico_barras <-
    puntos |>
    ggplot(aes(x = fct_infreq(t2_uso), fill = fct_infreq(t2_uso))) +
    geom_bar(
    aes(
      text = paste0(
        "Cantidad de puntos: ", after_stat(count)
      )     
    )
    ) +
    scale_fill_manual(values = COLORES_USO) +
    xlab("Uso") +
    ylab("Cantidad de puntos") +
    labs(fill = "Uso") +
    theme_clean() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.title = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank()
  )

# Gráfico de barras plotly
ggplotly(grafico_barras, tooltip = "text") |> 
  config(locale = 'es') |>
  layout(
    showlegend = TRUE,
    legend = list(
      title = list(text = ''),
      bordercolor = 'rgba(0,0,0,0)' # Establece el color del borde a transparente
    )
  ) 

Código para el despliegue del gráfico de pastel

# Total de puntos
conteo <- puntos |>
    count(t2_uso)

# Asegurarse de que los colores se asignen en el orden correcto
colores_para_plotly <- COLORES_USO[as.character(conteo$t2_uso)]

# Gráfico de pastel plotly
plot_ly(
    conteo, 
    labels = ~ t2_uso, 
    values = ~ n,
    marker = list(colors = colores_para_plotly)
) |>
add_pie()

Análisis con el paquete FIESTA

Cobertura de la tierra

Porcentaje de coberturas por parcela en t1

Se utiliza la función datPBpnt2pct para transponer los datos de puntos a porcentajes en cada parcela de monitoreo.

# Transposición de puntos a porcentajes de coberturas por parcela en t1
porcentaje_cobertura_parcela_t1 <- 
    datPBpnt2pct(puntos, uniqueid="plot_id", tvar="t1_cobertura")

Código para el despliegue de la tabla

# Despliegue de la transposición de puntos a porcentajes de coberturas por parcela en t1
porcentaje_cobertura_parcela_t1 |>
  datatable(
    caption = "Porcentaje de coberturas por parcela en t1",
    rownames = FALSE,
    options = list(
      pageLength = 5,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Porcentaje de coberturas por parcela en t2

# Transposición de puntos a porcentajes de coberturas por parcela en t2
porcentaje_cobertura_parcela_t2 <- 
    datPBpnt2pct(puntos, uniqueid="plot_id", tvar="t2_cobertura")

Código para el despliegue de la tabla

# Despliegue de la transposición de puntos a porcentajes de coberturas por parcela en t2
porcentaje_cobertura_parcela_t2 |>
  datatable(
    caption = "Porcentaje de coberturas por parcela en t2",
    rownames = FALSE,
    options = list(
      pageLength = 5,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Porcentaje total de coberturas en t1

La función modPBpop() genera datos de poblaciones.

# Generación de datos de poblaciones
PBpopdat <- modPBpop(pnt = puntos,
    pltassgnid = "plot_id",
    pntid = "sample_id")

La función modPB() genera las estimaciones.

# Estimación de la distribución de coberturas en t1
LCt1 <- modPB(PBpopdat = PBpopdat, rowvar = "t1_cobertura")

results.LCt1 <- LCt1$est

Código para el despliegue de la tabla

# Despliegue de la distribución de coberturas en t1
results.LCt1 |>
  datatable(
    caption = "Distribución de coberturas en t1",
    rownames = FALSE,
    options = list(
      pageLength = 5,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Porcentaje total de coberturas en t2

# Estimación de la distribución de coberturas en t2
LCt2 <- modPB(PBpopdat = PBpopdat, rowvar = "t2_cobertura")

results.LCt2 <- LCt2$est

Código para el despliegue de la tabla

# Despliegue de la distribución de coberturas en t2
results.LCt2 |>
  datatable(
    caption = "Distribución de coberturas en t2",
    rownames = FALSE,
    options = list(
      pageLength = 5,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Cambios en las áreas de coberturas de t1 a t2

PBpoparea <- 
    modPBpop(
      pnt = puntos,
      pltassgnid = "plot_id",
      pntid = "sample_id",
      unitarea = AREA_ESTUDIO
  )

# Cobertura en t1 vs cobertura en t2
coberT1vT2 <- 
    modPB(
        PBpopdat = PBpoparea,
      rowvar = "t1_cobertura",
      colvar = "t2_cobertura",
        tabtype = "AREA"
    )

Código para el despliegue de la tabla

# Despliegue de datos en una tabla
coberT1vT2$est |>
  datatable(
    caption = "Cobertura en t1 vs cobertura en t2",
    rownames = FALSE,
    options = list(
      pageLength = 6,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Flujo de cambios de cobertura de la tierra de t1 a t2

Código para el despliegue del gráfico de Sankey

puntos_sankey <-
    puntos |>
    select(sample_id, t1_cobertura, t2_cobertura)

# write_csv(puntos_sankey, "puntos_sankey.csv")
# puntos_sankey <- read_csv("puntos_sankey.csv")

# El prefijo es para diferenciar t1 y t2 en el gráfico de Sankey
puntos_sankey <- 
    puntos_sankey |>
  mutate(
    t1_cobertura = paste0("T1-", t1_cobertura),
    t2_cobertura = paste0("T2-", t2_cobertura)
  )

label <- c("T1-1000-Vegetacion", "T1-2000-Sin vegetacion", "T1-3000-Agua", 
                     "T1-4000 Nubes y sombras", "T1-5000-Sin informacion", "T1-Nulo",
                     "T2-1000-Vegetacion", "T2-2000-Sin vegetacion", "T2-3000-Agua", 
                     "T2-4000 Nubes y sombras", "T2-5000-Sin informacion", "T2-Nulo"
         )

color <- c(COLOR_VEGETACION, COLOR_SIN_VEGETACION, COLOR_AGUA, 
           COLOR_NUBES_Y_SOMBRAS, COLOR_SIN_INFORMACION, COLOR_NULO,
           COLOR_VEGETACION, COLOR_SIN_VEGETACION, COLOR_AGUA, 
           COLOR_NUBES_Y_SOMBRAS, COLOR_SIN_INFORMACION, COLOR_NULO
         )

# Agregar columnas de códigos para usar en el gráfico de Sankey
puntos_sankey <- puntos_sankey |>
  mutate(
    t1_codigo = case_when(
      t1_cobertura == 'T1-1000-Vegetacion' ~ 0,
      t1_cobertura == 'T1-2000-Sin vegetacion' ~ 1,
      t1_cobertura == 'T1-3000-Agua' ~ 2,
      t1_cobertura == 'T1-4000 Nubes y sombras' ~ 3,
      t1_cobertura == 'T1-5000-Sin informacion' ~ 4,
      t1_cobertura == 'T1-Nulo' ~ 5
    ),
    t2_codigo = case_when(
      t2_cobertura == 'T2-1000-Vegetacion' ~ 6,
      t2_cobertura == 'T2-2000-Sin vegetacion' ~ 7,
      t2_cobertura == 'T2-3000-Agua' ~ 8,
      t2_cobertura == 'T2-4000 Nubes y sombras' ~ 9,
      t2_cobertura == 'T2-5000-Sin informacion' ~ 10,
      t2_cobertura == 'T2-Nulo' ~ 11
    )
  )

# Contar las combinaciones de t1_codigo y t2_codigo
conteo_combinaciones <- 
    puntos_sankey |>
  count(t1_codigo, t2_codigo)

# Crear las listas separadas
source <- conteo_combinaciones$t1_codigo
target <- conteo_combinaciones$t2_codigo
value <- conteo_combinaciones$n

# Gráfico de Sankey
grafico_sankey <- 
    plot_ly(
    type = "sankey",
    orientation = "h",
    node = list(
      label = label,
      color = color,
      pad = 15,
      thickness = 20,
      line = list(
        color = "black",
        width = 0.5
      ),
      font = list(
        size = 12,
        color = "#000000",
        weight = "bold"
      )
    ),
    link = list(
      source = source,
      target = target,
      value =  value
    )
  )

grafico_sankey <- 
    grafico_sankey |> 
    layout(
    font = list(
      size = 12,
      color = "#000000",
      weight = "bold"
    )  
  )

grafico_sankey

Uso de la tierra

Porcentaje de usos por parcela en t1

# Transposición de puntos a porcentajes de usos por parcela en t1
porcentaje_uso_parcela_t1 <- 
    datPBpnt2pct(puntos, uniqueid="plot_id", tvar="t1_uso")

Código para el despliegue de la tabla

# Despliegue de la transposición de puntos a porcentajes de uso por parcela en t1
porcentaje_uso_parcela_t1 |>
  datatable(
    caption = "Porcentaje de usos por parcela en t1",
    rownames = FALSE,
    options = list(
      pageLength = 5,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Porcentaje de usos por parcela en t2

# Transposición de puntos a porcentajes de uso por parcela en t2
porcentaje_uso_parcela_t2 <- 
    datPBpnt2pct(puntos, uniqueid="plot_id", tvar="t2_uso")

Código para el despliegue de la tabla

# Despliegue de la transposición de puntos a porcentajes de usos por parcela en t2
porcentaje_uso_parcela_t2 |>
  datatable(
    caption = "Porcentaje de usos por parcela en t2",
    rownames = FALSE,
    options = list(
      pageLength = 5,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Porcentaje total de usos en t1

# Generación de datos de poblaciones
PBpopdat <- modPBpop(pnt = puntos,
    pltassgnid = "plot_id",
    pntid = "sample_id")

# Estimación de la distribución de usos en t1
LCt1 <- modPB(PBpopdat = PBpopdat, rowvar = "t1_uso")

results.LCt1 <- LCt1$est

Código para el despliegue de la tabla

# Despliegue de la distribución de usos en t1
results.LCt1 |>
  datatable(
    caption = "Distribución de usos en t1",
    rownames = FALSE,
    options = list(
      pageLength = 5,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Porcentaje total de usos en t2

# Estimación de la distribución de usos en t2
LCt2 <- modPB(PBpopdat = PBpopdat, rowvar = "t2_uso")

results.LCt2 <- LCt2$est

Código para el despliegue de la tabla

# Despliegue de la distribución de usos en t2
results.LCt2 |>
  datatable(
    caption = "Distribución de usos en t2",
    rownames = FALSE,
    options = list(
      pageLength = 5,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Cambios en las áreas de usos de t1 a t2

PBpoparea <- 
    modPBpop(
      pnt = puntos,
      pltassgnid = "plot_id",
      pntid = "sample_id",
      unitarea = AREA_ESTUDIO
  )

# Uso en t1 vs uso en t2
usoT1vT2 <- 
    modPB(
        PBpopdat = PBpoparea,
      rowvar = "t1_uso",
      colvar = "t2_uso",
        tabtype = "AREA"
    )

Código para el despliegue de la tabla

# Despliegue de datos en una tabla
usoT1vT2$est |>
  datatable(
    caption = "Uso en t1 vs uso en t2",
    rownames = FALSE,
    options = list(
      pageLength = 6,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Flujo de cambios de uso de la tierra de t1 a t2

Código para el despliegue del gráfico de Sankey

puntos_sankey <-
    puntos |>
    select(sample_id, t1_uso, t2_uso)

# write_csv(puntos_sankey, "puntos_sankey.csv")
# puntos_sankey <- read_csv("puntos_sankey.csv")

# El prefijo es para diferenciar t1 y t2 en el gráfico de Sankey
puntos_sankey <- 
    puntos_sankey |>
  mutate(
    t1_uso = paste0("T1-", t1_uso),
    t2_uso = paste0("T2-", t2_uso)
  )

label <- c("T1-1000-MC bosques", "T1-2000-Agricultura", "T1-3000-Ganad y past", 
                     "T1-4000-Zonas humed", "T1-5000-Infraest", "T1-6000-Otras tierras",
                     "T1-7000 No clasif", "T1-Nulo",
                     "T2-1000-MC bosques", "T2-2000-Agricultura", "T2-3000-Ganad y past", 
                     "T2-4000-Zonas humed", "T2-5000-Infraest", "T2-6000-Otras tierras",
                     "T2-7000 No clasif", "T2-Nulo"
                    )

color <- c(COLOR_MC_BOSQUES, COLOR_AGRICULTURA, COLOR_GANADERIA_Y_PASTOS, 
           COLOR_ZONAS_HUMEDALES, COLOR_INFRAESTRUCTURA, COLOR_OTRAS_TIERRAS,
           COLOR_NO_CLASIFICACION, COLOR_NULO,
           COLOR_MC_BOSQUES, COLOR_AGRICULTURA, COLOR_GANADERIA_Y_PASTOS, 
           COLOR_ZONAS_HUMEDALES, COLOR_INFRAESTRUCTURA, COLOR_OTRAS_TIERRAS,
           COLOR_NO_CLASIFICACION, COLOR_NULO
                    )

# Agregar columnas de códigos para usar en el gráfico de Sankey
puntos_sankey <- puntos_sankey |>
  mutate(
    t1_codigo = case_when(
      t1_uso == 'T1-1000-MC bosques' ~ 0,
      t1_uso == 'T1-2000-Agricultura' ~ 1,
      t1_uso == 'T1-3000-Ganad y past' ~ 2,
      t1_uso == 'T1-4000-Zonas humed' ~ 3,
      t1_uso == 'T1-5000-Infraest' ~ 4,
      t1_uso == 'T1-6000-Otras tierras' ~ 5,
      t1_uso == 'T1-7000 No clasif' ~ 6,
      t1_uso == 'T1-Nulo' ~ 7
    ),
    t2_codigo = case_when(
      t2_uso == 'T2-1000-MC bosques' ~ 8,
      t2_uso == 'T2-2000-Agricultura' ~ 9,
      t2_uso == 'T2-3000-Ganad y past' ~ 10,
      t2_uso == 'T2-4000-Zonas humed' ~ 11,
      t2_uso == 'T2-5000-Infraest' ~ 12,
      t2_uso == 'T2-6000-Otras tierras' ~ 13,
      t2_uso == 'T2-7000 No clasif' ~ 14,
      t2_uso == 'T2-Nulo' ~ 15
    )
  )

# Contar las combinaciones de t1_codigo y t2_codigo
conteo_combinaciones <- 
    puntos_sankey |>
  count(t1_codigo, t2_codigo)

# Crear las listas separadas
source <- conteo_combinaciones$t1_codigo
target <- conteo_combinaciones$t2_codigo
value <- conteo_combinaciones$n

# Gráfico de Sankey
grafico_sankey <- 
    plot_ly(
    type = "sankey",
    orientation = "h",
    node = list(
      label = label,
      color = color,
      pad = 15,
      thickness = 20,
      line = list(
        color = "black",
        width = 0.5
      ),
      font = list(
        size = 12,
        color = "#000000",
        weight = "bold"
      )
    ),
    link = list(
      source = source,
      target = target,
      value =  value
    )
  )

grafico_sankey <- 
    grafico_sankey |> 
    layout(
    font = list(
      size = 12,
      color = "#000000",
      weight = "bold"
    )
  )

grafico_sankey